#include "../Random.hpp"
#include "../Kernel.hpp"
//Implementation of Robert's two sided truncated Gaussian (1995)
class Rtnorm
{
	public:
	Rtnorm(void)
	{
		boost::random::uniform_01<> U;
		runif=new RandomG::Random<boost::random::uniform_01<> >(U);
	}

	double operator()(double m,double mm,double mp,double s)
	{
		double res=0;
		double i=0;
		mm=mm-m;
		mp=mp-m;
		while(i==0)
		{
			double z=(mp-mm)*(*runif)()+mm;
			double r=rho(z,mm,mp);	
			if(u<=r)
			{
				i=1;
				res=z*s+m;				
			}
		}
	}
	double rho(double z, double mm, double mp)
	{
		double res=0;
		if(0>mp)
		{
			res=exp((pow(mp,2)-pow(z,2))/2);
		}else if(mm>0){
			res=exp((pow(mm,2)-pow(z,2))/2);
		}else{
			res=exp(-pow(z,2)/2);
		}
		return res;
	}

	private:
	RandomG::Random<boost::random::uniform_01<> > *runif;
}



